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The existence of universal scaling in the vicinity of the jamming transition of sheared 
granular materials is predicted by a phenomenology. The critical exponents are explicitly 
determined, which are independent of the spatial dimension. The validity of the theory is 
verified by the molecular dynamics simulation. 

Jamming is an athermal phase transition between the sohd-Hke jammed phase 
and the hquid-hke unjammed phase of granular assembhes. Above the critical den- 
sity, which is referred to the point J, the assemblies obtain the rigidity and the 
dynamic yield stress, while assemblies behave like dense liquids below the point J. 
Liu and Nagel^^ indicated that the jamming transition is a key concept of glassy 
materials. Since then many aspects of similarities between the conventional glass 
transition and the jamming transition have been investigated Indeed, there are 
many examples where granular materials are used in order to investigate dynamical 
heterogeneity in glassy materials.'^ '•^•'''^''^ On the other hand, we still do not have 
an unified view in describing glassy materials because we cannot use the conventional 
theoretical tool for the glass transition such as the mode-coupling theory.l^ 

The jamming is a continuous transition that bulk and shear moduli become 
nonzero, and there are scaling laws in the vicinity of the point J similar to the 
cases of conventional critical phenomena-'^!'® Olsson and TeiteP^ and HatanJIS 
further demonstrated the existence of beautiful scalings near the point J. Therefore 
we can expect the existence of a simple theory in describing the jamming transition. 
However, we still do not have any theory to determine the critical exponents of the 
jamming transition. 

In this letter, we predict the critical exponents for jamming transition of sheared 
granular materials based on a phenomenology. First, we introduce the system we 
consider and the scaling laws, some of which were introduced in Ref. 11). Second, 
based on a phenomenological theory, we decide the critical exponents for the jamming 
transition. Finally, we verify the theoretical prediction from our simulation. The 
most surprising in our finding is that the critical exponents are independent of the 
spatial dimension D, but depend on the type of particle interaction. 

Let us consider a dense sheared and frictionless granular system in which uniform 
shear flow is stable. The system consists of N spherical grains in D dimensions. An 
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important parameter to characterize the system is the volume fraction (p. In contrast 
to granular gases, the contact force plays crucially important roles in jamming tran- 
sition, where the normal contact force is repulsive one characterized by k5^. Here, 
k is the stiffness constant, and 612 = r — (di + o"2)/2 with the relative distance r be- 
tween the contacting particles of the diameters cJi and a2- We believe that Hertzian 
contact law A = 3/2 is appropriate for three dimensional grains, but we often use 
a simpler linear spring model with A = 1. In this letter, we omit any tangential 
contact force between grains. Thus, granular particles are frictionless, which can 
simplify the argument. 

As introduced in Ref. [TT]) . this system is expected to exhibit the scalings for 
the granular temperature T and the shear stress S in the sheared plane near the 
jamming transition as 

T = Ar,n\<^m{tn^), (1) 

where 7 is the shear rate, (p = (f) — is the excess volume fraction from the critical 
fraction cpj at the point J, At^d and As^d are respectively amplitudes of the granular 
temperature and the shear stress. The scaling functions T^{x) and S^{x) above the 
point J respectively differ from T-{x) and S-{x) below the point J. It should be 
noted that At^d, As^Di to and do not depend on 7 and c^, but depend only on 
D. 

We also introduce a characteristic time scale of sheared granular assemblies as 

P> 

where n is the number density of grains. This oj is reduced to the collision frequency 
in the unjammed phase in a steady state achieved by the balance between the viscous 
heating and the collisional energy loss. In contrast to the assumption by Hatano et 
al.^ to also satisfies the scaling form 

u; = A^M<l>rw^[^Dj^). (4) 

Thus, there are six critical exponents x^, x^, y^, y-y, z<p, and in eqs. ([T|), ([2]) 
and dH). We note that the normal stress P also satisfies a similar scaling relation 
P = Ap^£)|(^|^*7^-|-(p£)7/|^P t),--* but we omit the details of the arguments on 
y' ^ and y' ^ in this letter. We will discuss them elsewhere. 

Below the point J, Bagnold's scaling should be held.'^''^''^ Thus, the scaling 
functions in the unjammed branch satisfy 

lim7L(x) = lim5_(x) = x^, limW_(x) = x. (5) 

x— >0 a;— >0 x— >0 

On the other hand, the jammed branch is characterized by the dynamic yield stress 
and the freezing of motion. Then, the scaling functions in the jammed branch satisfy 

lim = lim >V+(x) = 1, lim '7+(x) = x. (6) 

x— >0 x^O x— >0 
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To obtain the last equation in eq. ([6]), we have used eq. ([3]) and the other two 
equations in eq. ([6]). Since the scahng functions are independent of ^ at the point 
J, we obtain 

hm T±{x) oc x^^, hm S±{x) tx x^^, hm W±{x) tx x^^ . (7) 

X— »oo X— ►oo z— >oo 

Now, let us determine the six critical exponents. First, we note that there are 
three trivial relations among the exponents. From eqs. and ^ we obtain 

Zcp = y^p - x^{l - x~^). (8) 

Similarly, from eqs.([T|)-(|l]) with eq. ([5]) or ([7]) we respectively obtain 

z<p{l - z-^) = y^{l - 2y-i) - x^(l - 2x~^), (9) 
= — x^ + 1. (10) 

Thus, we further need three relations to determine the exponents. 

In order to introduce the other relations, we consider the pressure P in the limit 
7 — > 0. Let us consider Cauchy's stress in the jammed phase in which the pressure P 
is given by P = Ylii>j ■> where V is the volume of the system, rij and Fij are 

respectively the distance and the force between i and j particles. This expression may 
be approximated by P ~ Z{^)rc{^)Fc{^) in the zero shear limit, where Z{$) is the 
average coordination number, rc{^) and Fc{^) are respectively the average distance 
between contacting grains and the average force acting on the contact point. It is 
obvious that Z{<P) and rc(<P) can be replaced by Z{Q) and rc(0) = a in the vicinity 
of the jamming point, where a is the average diameter of the particles. Indeed, 
O'Hern et a\.¥^ verified Z{^) — Z{{)) oc ^^/-^ for three dimensional cases. Thus, the 
most important term is the mean contact force Fc{^) where S{'P) is the 

average length of compression. Now, let us compress the system at the critical point 
(f)j into ^ = (j) — > by an affine transformation. Since all the characteristic 
lengths are scaled by the system size, we may assume the approximate relation 
rc(^) = / a. From the relation 5{(P) = ^^(O) — rci^P), approximately 
satisfies 5{<P) ~ {a/D4>j)0 ~ ^ in the vicinity of ^ = 0. Thus, we conclude P ~ 'P^. 
This relation has also been verified in Ref . [9]) . 

On the other hand, it is well-known that there is Coulomb's frictional law in gran- 
ular systems in which S/P is a constant. Indeed, Hatand^® simulated the sheared 
granular system under a constant pressure P and demonstrated that the ratio sat- 
isfies limAf^oS{^,P)/P = Sy{P)/P = Mo, where 5y(P) = lim^_o5(7,P) and the 
constant Mq is independent of the pressure P. The excess volume fraction ?P(7, P) in 
this system is a function of 7 and P, but we also can express the pressure as P(7, P). 
Since Mq = lim^^o 5'(7, P(7, ^))/P(7, ^) = 5^(^(0, cP))/P(0,<?) is independent of 
P(0,^), Mo should be independent of (p. Thus, we can conclude that ^ dependence 
of S is the same as that of P in the limit 7 — > 0. From this result, we obtain y,p = A 
and 

2/^ = 1 for linear spring model. (11) 

This result also implies y$ = y' ^ which is consistent with the numerical observa- 
tion.EJ 
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The next relation is related to the density of state. Wyart et al.^"^' demonstrated 
the followings for unsheared assemblies of elastic soft spheres, (i) The jamming is 
related to the appearance of the soft modes in the density of state, (ii) There is 
a plateau in the density of state in the vicinity of jamming transition, (iii) The 
cutoff frequency w* of the plateau is proportional to ^/P. From the argument in the 
previous paragraph, the applied pressure satisfies the relation P oc 'P^. When we 
assume that the characteristic frequency uj in the limit 7 — > can be scaled by the 
cutoff frequency uj*, we may conclude lo ~ |(^|^/^ for the linear spring model. Thus, 
we obtain = A/2 or 

Z4, = 1/2 for linear spring model. (12) 

Finally, we consider the characteristic frequency uj in the unjammed phase {<P < 
0). In this phase, the characteristic frequency uj is estimated as a; ~ y^T/'m/l(<P), 
where 1{P) is the mean free path. Note that 1{P) may be evaluated as ((j/D(/>j)|^| 
in the vicinity of the point J, using the parallel argument to S{P) for > 0. 
From the scalings in Bagnold's regime ([5]), we obtain u! ~ [^[^^(i"^^ )j and T ~ 
^^^X4,(i-2x-, )^2^ Substituting these relations to a; ~ ^yT/m/l{P), we obtain 

z<f (1 - z^i) - ^x^(l - 2x7!) = -1. (13) 

From the above six relations ([8])- ()13p we finally determine the six critical expo- 
nents 
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y-y = -^, 2^ ^'^ " 5 ^^^^ 

for the linear spring model. The exponents in Hertzian model are, of course, different. 
In general situation for A, eqs. ()14p are replaced by 

2A + 4 

x<p = 2 + A, x^= y<P = A, 

2A A A 

^^ = ZT4' '''=r '^ = aTI- ^^^^ 

We should note that the exponents are independent of the spatial dimension. This 
is not surprising because our phenomenology to derive eqs. (fTTj) and (fT2]) is inde- 
pendent of the spatial dimension.lS'l^i "\Ye should stress an interesting feature of 
jamming transition that the exponents strongly depend on the interaction model 
among particles. This property is contrast to that in the conventional critical phe- 
nomena. Thus, we should be careful to use the idea of the universality in describing 
the jamming transition. 

From now on, let us verify our theoretical results based on the molecular dynam- 
ics simulation. In our simulation, the system consists of N spherical grains in 2, 3, 4 
dimensions. We adopt the linear spring model ( Z\ = 1 ) for simplicity. We also 
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introduce dissipative force —rjdv, where 5v represents the relative velocity between 
the contacting particles. Each grain has an identical mass m. In order to realize 
an uniform velocity gradient 7 in y direction and macroscopic velocity only in the x 
direction, we adopt the Lees-Edwards boundary conditions. The particle diameters 
are 0.7c70) O.Scto, 0.9(To and do each of which is assigned to N/A particles. 

In our simulation m, do and rj are set to be unity, and all quantities are converted 
to dimensionless forms, where the unit of time scale is m/rj. We use the spring 
constant k = 1.0. For the system near the critical density, such as 4> = 0.8428 for 
D = 2, (j) = 0.643 and 0.6443 for D = 3, we use N = 4000 in order to remove finite 
size effects. For other systems, we use N = 2000. 

The scaling plots of our simulation based on the exponents (jl4p are shown in 
Fig. [TJ We should stress that these scaling plots contain the data in D = 2, 3 and 
4. The volume fraction at the point J is estimated as (pj = 0.84285 for D = 2, 
4>j = 0.64455 for D = 3 01 = 0.4615 for D = 4. We examine the shear rate 7 is 
in the range between 5 x 10~^ and 5 x 10"^ for D = 2,3 and between 5 x 10~® and 
5 X 10~^ for D = 4. The amplitudes and the adjustable parameters are obtained 
as {tn, At,D, SD, As,D,WD, A^^d) = (0.0125,7.17,0.025,0.035,0.05,0.3) for D = 2, 
(0.01385, 2.527, 0.03, 0.04, 0.065, 0.65) for D = 3, (0.015, 1.6275, 0.03, 0.06, 0.06, 1) for 
D = 4. Since Fig {1] exhibits beautiful scaling laws, our phenomenology seems to be 
right. 

Figure 2(a) shows ^ dependence of the shear viscosity ^ = S/'y in Bagnold's 
regime. We should note that there is no consensus in previous studies on the shear 
viscosity. For example, Garcia-Rojo et alW^ reported fi ~ l/((/'c — 0)) where (pc is 
lower than while Losert et al.'^ observed the exponent larger than 1 from their 
experiment, and the exponent of divergence in Ref. [TOl) is also between 1 and 2. We 
also note that the viscosity is believed to diverge as |^^|~^ for colloidal suspensions 
However, our scaling theory predicts ^ oc \<p\y't>(^~'^/y-i) esc |^^|~^, and the viscosity 
diverges at the point J. Here, the scaling exponents for is independent of A because 
1/(1) and for an arbitrary A are determined from our theory as in eq. (|15p . As we 
can see in Fig. 2(a), the theoretical prediction is consistent with our numerical 
result. We also examine the possibility that the viscosity diverges at (pc < (f'J with 
fJ- ~ {'Pc — 4')~^ in the case of Garcia-Rojo et al.'^^ Actually we can fit the data 
of our two-dimensional simulation by ~ (0c — 0) ""^ with cp^ = 0.835 which is less 
than (f)j = 0.8428 for (p < (j)c, but the viscosity is still finite even for (p > (j)c (see Fig. 
2(b)). Thus, we can conclude that (i) the viscosity does not satisfy {(pc — (p)~^ but 
exhibits a consistent behavior with {cpj — (p)~'^ predicted by our phenomenology, and 
(ii) the critical behaviors are only characterized by the point J. 

We also verify the validity of w ~ |(^| in the jammed phase from our simulation 
in Fig. 3. The envelope line of our result seems to be consistent with the theoretical 
prediction. 

It should be noted that the existence of the plateau for \^\ ^ in Fig. 2(a) 
can be understood from the scaling relation ([2]). Indeed, ^ dependence of the shear 
stress disappears in the limit of large x = 'y/\^\y*/y'i with S-{x) — > x^'' as in eq.©. 
Thus, we obtain fx/j = S/j'^ ~ jy''~'^ ~ 7^®/^. This estimation might be consistent 
with the simulation in Fig. 2(a) in which the value of the plateau increases as the 
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Fig. 1. (Color online) (a) : Collapsed data of the shear rate dependence of the granular temperature 
T using the scaling law for D — 2,3 and 4. The dashed line, the dotted line and the solid line 
are proportional to 7, 7^ and 7^^ . The legends show the dimension D and the volume fraction 
(j) as {D,(j)). (b): Collapsed data of the shear rate dependence of the shear stress S using the 
scaling law for D — 2,3 and 4. The dotted line and the solid line are proportional to 7^ and 
7^^. (c): Collapsed data of the shear rate dependence of the cooling rate u using the scaling 
law for D = 2,3 and 4. The dotted line and the solid line are proportional to 7'^ and 7^^. 
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shear rate decreases. Similarly, we can expect the value of plateau of uj as , while 
this saturation cannot be verified from the simulation. 
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Fig. 2. (a) : ^/-y as a function of for D = 2 with 7 = 5 x 10"^ 5 x 10~®,5 x 10"'^ in the 
unjammed phase, (b) : /i/7 as a function of 4> for 7 = 5 x 10^^ in the unjammed phase, where 
the solid line is proportional to {(t>c — <t>)~^ with (f>c = 0.835. 




Fig. 3. UJ as a. function of $ for D = 2 with 7 = 5 x 10 ®, 5 x 10 ^, 5 x 10 * in the jammed phase. 

Now, let US discuss our results. First of all, the ratios between the exponents 
y^/y^, and z^/z^ obtained in eq. ([HD or satisfy 

a^^ = ^ = ^ = ^. (16) 

This is not surprising because the time scale is expected to be scaled by the shear 
rate. Thus, the ratio a in eqs. ([I]), ([1]), and @ should be common. In other words, 
the characteristic time scale r exhibits the critical slowing down as r ~ |^|~"- This 
property has already been indicated by Hatano.'^I' Once we accept the ansatz (fT6|) . 
eqs. dS]), dH), and ([10]) are degenerate, and reduce to 



X(j,-y(j, + z^ = a. 



(17) 
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Equation (I13p is also reduced to the simplified form 

^0 = f - 1- (18) 
From these equations and eqs. ([TT]) and (fT2]) with p^ . we obtain eq. (fT^ or eq. 

(HSI). 

Second, Hatano estimated the exponents = 2.5, = 1.3, = 1.2 and 
y-y = 0.57 from his three-dimensional simulation for the linear spring model^^^' which 
differ from our prediction (jl4p . In particular, if we use these values with eq. (jl8p . 
Zfj, is estimated as = 0.25, which is one half of our prediction z<p = 1/2. However, 
the estimation of the scaling exponents strongly depend on the choice of and the 
range of the shear rate 7. The value of and the range of 7 in Ref. [TT]) are larger 
than ours. If we adopt Hatano's (pj and the range of 7, our numerical data can be 
scaled by Hatano's scaling. Although his scaling can be used in the wide range of 7, 
the deviation from his scaling can be detected in the small 7 region (7 < 10~^). It 
is obvious that we should use smaller 7 as possible as we can to extract the critical 
properties. This suggests that our exponents are more appropriate than Hatano's 
exponents in characterizing the jamming transition. The difficulty in determination 
of the exponents from the simulation also supports the significance of our theory to 
determine the scaling laws. 

Third, the exponents obviously depend on the model of interaction between 
particles as predicted in eq. (fl5]) . Our preliminary simulation suggests that the 
numerical exponents for Hertzian contact model is consistent with the prediction 
of ()15p . The numerical results on A dependence of the exponents will be reported 
elsewhere. 

Fourth, our results should be modified when we analyze the model in the zero 
temperature limit of Langevin thermostat. This situation corresponds to that in 
Ref. nop . In this case, we should replace Bagnold's law in unjammed phase by 
Newtonian law oc 7. As a result, all the scaling exponents have different values. 
We will discuss the results of this situation elsewhere. 

Finally, we comment on the relation between our results and the previous studies 
on dynamical heterogeneity in glassy materials. The dynamical heterogeneity in 
glassy materials is characterized by the large fluctuations of four point correlation 
function, in which the result strongly depends on the spatial dimension. On the other 
hand, our theory and numerical simulation suggest that the critical fluctuation is not 
important and our phenomenology works well. Since the quantities we analyzed in 
this letter are not directly related to the four-point correlation functions, there is no 
distinct contradiction between them. To study the roles of critical fluctuations and 
dynamical heterogeneity we may need a more sophisticated theory. This will be our 
future task. 

In conclusion, we develop the phenomenological theory in describing the jam- 
ming transition. We determine the critical exponents which are independent of the 
spatial dimension. The validity of our theory has been verified by the molecular 
dynamics simulation. 
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